Age : Age of the patient
Sex : Sex of the patient
exang: exercise induced angina (1 = yes; 0 = no)
ca: number of major vessels visible in fluoroscopy (0-3)
cp : Chest Pain type chest pain type Value 1: typical angina Value 2: atypical angina Value 3: non-anginal pain Value 4: asymptomatic
trtbps : resting blood pressure (in mm Hg)
chol : cholestoral in mg/dl fetched via BMI sensor
fbs : (fasting blood sugar > 120 mg/dl) (1 = true; 0 = false)
rest_ecg : resting electrocardiographic results Value 0: normal Value 1: having ST-T wave abnormality (T wave inversions and/or ST elevation or depression of > 0.05 mV) Value 2: showing probable or definite left ventricular hypertrophy by Estes' criteria
thalach : maximum heart rate achieved
thal : Thalium Stress Test result
slope : the slope of the peak exercise ST segment (2 = upsloping; 1 = flat; 0 = downsloping)
oldpeak : ST depression induced by exercise relative to rest
Then, calculate the what-if explanations of these predictions using Ceteris Paribus profiles (also called What-if plots), e.g. in Python: AIX360, Alibi dalex, PDPbox; in R: pdp, DALEX, ALEplot. *implement CP yourself for a potential bonus point


Interesting thing here is that for both observations the blood pressure that would minimize the rist of getting the heart attack is higher then the recommended (around 120).
It is especially visible for the second observation. I think model could kind of underestimate the importance of the resting blood pressure, because it is strongly correleted to the other variables.
Find two observations in the data set, such that they have different CP profiles. For example, model predictions are increasing with age for one observation and decreasing with age for another one. NOTE that you will need to/ have a model with interactions to observe such differences.
Such situation is happening for the observations I have mentiond above.
First observations is a woman, that is pretty healthy, so her rist of getting heart attack increases with age as we would expect.
In the second observations we have a rather unhealthy man, who at a not very old age have some bad examination results such as very elevated cholesterols. Such thing is pretty rare in younger people, so it might mean that his health problems are very serious, while at older age it could be just the sign of normal decrease in healthiness.
Compare CP, which is a local explanation, with PDP, which is a global explanation. *implement PDP yourself for a potential bonus point

Here we can see the global explanations. I think they are pretty similar to what we have seen in the the local explanations. We can see that differences betwwen maximum and minimum are pretty small on both plots.
The resting blood pressure plot is very suprising, since all the medical research says that high resting blood pressure is a bad sign. There are a few potential explanated to why this plot might look like this:
Compare PDP between between at least two different models.

We can see that the plots look different then before. It is because logistic regression is a linear model.
The resting heart pressure have similar tendency as before (heart attack risk is decreasing with increasment in resting heart pressure).
For age the profile is strongly differetnt. Here we just get an increasing plot, while before we had much more complex relation.
In the report I am still using the implementations from dalex, because they look better. (The plots are labeld and there is nice grid in a back. It is pretty easy to add this in matplotlib, but I wanted to keep code very simple and straightforward.)
import matplotlib.pyplot as plt
def CP(model, observation, variables, minVals, maxVals):
for i, variable in enumerate(variables):
x = np.linspace(minVals[i],maxVals[i], num=75)
observations = []
predictions = []
for j in x:
o = observation
o[[variable]] = j
predictions.append(model.predict_proba(o)[:,1])
predictions.append(model.predict_proba(o)[:,1])
x = np.append(x, observation[variable])
plt.plot(x, predictions)
plt.show()
import copy
def PDP(model, X, variable, minVal, maxVal):
predictions = []
dataset = copy.deepcopy(X)
dataset['age'] = 5
x = np.linspace(minVal, maxVal, num=75)
for i, point in enumerate(x):
predictions.append([])
dataset[variable] = point
predictions[i].append(model.predict_proba(dataset)[:,1])
predictions = np.array(predictions)
predictions = predictions.mean(axis=2)
plt.plot(x, predictions)
plt.show()

import pandas as pd
import sklearn
from sklearn import ensemble
import dalex as dx
import lime
dataset = pd.read_csv('heart.csv')
dataset = pd.get_dummies(dataset)
dataset
features = dataset.drop(columns='output')
#fixing typo in data
features['thalach']=features['thalachh']
features = features.drop(columns='thalachh')
features['slope']=features['slp']
features = features.drop(columns='slp')
features['ca']=features['caa']
features = features.drop(columns='caa')
features = pd.get_dummies(features, columns=['cp', 'thall'])
features
X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(features, dataset['output'], test_size=0.3, random_state=0)
X_train
forest = sklearn.ensemble.RandomForestClassifier()
forest.fit(X=X_train,y=y_train)
print(f'Accuracy: {sklearn.metrics.accuracy_score(y_test,forest.predict(X_test))}')
print(f'Recall: {sklearn.metrics.recall_score(y_test,forest.predict(X_test))}')
print(f'Precision: {sklearn.metrics.precision_score(y_test,forest.predict(X_test))}')
forest_accuracy = sklearn.metrics.accuracy_score(y_test,forest.predict(X_test))
forest_recall = sklearn.metrics.recall_score(y_test,forest.predict(X_test))
forest_precision = sklearn.metrics.precision_score(y_test,forest.predict(X_test))
print('\nResults on train dataset:')
print(f'Accuracy: {sklearn.metrics.accuracy_score(y_train,forest.predict(X_train))}')
print(f'Recall: {sklearn.metrics.recall_score(y_train,forest.predict(X_train))}')
print(f'Precision: {sklearn.metrics.precision_score(y_train,forest.predict(X_train))}')
obs = list(range(2))
import numpy as np
predict = lambda m, d: m.predict_proba(d)[:,1]
for i in obs:
obs1 = X_test.iloc[obs[i]].to_numpy().reshape(1,-1)
print(obs1)
print(predict(forest, obs1))
X = X_test
y = y_test
explainer = dx.Explainer(forest, X_test, y_test, predict_function=predict, label="forest")
explainer.model_performance(cutoff=y.mean())
explainer.model_parts().result
CP(forest, X.iloc[[0]], ['age', 'trtbps'], [min(X['age']), min(X['trtbps'])], [max(X['age']), max(X['trtbps'])])
PDP(forest, X=X, variable='age', minVal=min(X['age']), maxVal=max(X['age']))
cp = explainer.predict_profile(new_observation=X.iloc[[0]])
cp.plot(variables=["age", "trtbps"])
cp = explainer.predict_profile(new_observation=X.iloc[[11]])
cp.plot(variables=["age", "trtbps"])
dpd = explainer.model_profile()
dpd.plot(variables=["age", "trtbps"])
import sklearn.linear_model
model = sklearn.linear_model.LogisticRegression(max_iter=500)
model.fit(X=X_train,y=y_train)
print(f'Accuracy: {sklearn.metrics.accuracy_score(y_test,model.predict(X_test))}')
print(f'Recall: {sklearn.metrics.recall_score(y_test,model.predict(X_test))}')
print(f'Precision: {sklearn.metrics.precision_score(y_test,model.predict(X_test))}')
model_accuracy = sklearn.metrics.accuracy_score(y_test,model.predict(X_test))
model_recall = sklearn.metrics.recall_score(y_test,model.predict(X_test))
model_precision = sklearn.metrics.precision_score(y_test,model.predict(X_test))
print('\nResults on train dataset:')
print(f'Accuracy: {sklearn.metrics.accuracy_score(y_train,model.predict(X_train))}')
print(f'Recall: {sklearn.metrics.recall_score(y_train,model.predict(X_train))}')
print(f'Precision: {sklearn.metrics.precision_score(y_train,model.predict(X_train))}')
X = X_test
y = y_test
predict = lambda m, d: m.predict_proba(d)[:,1]
explainer = dx.Explainer(model, X_test, y_test, predict_function=predict, label="lr")
print(explainer.model_performance(cutoff=y.mean()))
cp = explainer.model_profile()
cp.plot(variables=["age", "trtbps"])